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Abstract 

Finite element method is one of powerful numerical methods to solve PDE. Usu- 
ally, if a finite element solution to a Poisson equation based on a triangulation of the 
underlying domain is not accurate enough, one will discard the solution and then refine 
the triangulation uniformly and compute a new finite element solution over the refined 
triangulation. It is wasteful to discard the original finite element solution. We propose 
a prewavelet method to save the original solution by adding a prewavelet subsolution to 
obtain the refined level finite element solution. To increase the accuracy of numerical 
solution to Poisson equations, we can keep adding prewavelet subsolutions. 

Our prewavelets are orthogonal in the H 1 norm and they are compactly supported 
except for one globally supported basis function in a rectangular domain. We have 
implemented these prewavelet basis functions in MATLAB and used them for numerical 
solution of Poisson equation with Dirichlet boundary conditions. Numerical simulation 
demonstrates that our prewavelet solution is much more efficient than the standard 
finite element method. 



1 Introduction 

Finite element method is one of powerful numerical methods to solve PDE. Usually, if a finite 
element solution to a Poisson equation based on one level triangulation of the underlying 
domain is not accurate enough, one will discard the solution and then refine the triangulation 
and compute a new finite element solution at the refined level. It is wasteful to throw the 
original finite element solution away. In order to save the original solution and get the more 
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accurate new solution, we have to add H 1 orthogonal subsolution. That is, let V/, be a 
finite element space over a triangulation A/, and Vh/2 be the finite element space over the 
refined triangulation. Since Vh C Vh/2, let Wh = Vh/2 Q Vh under H 1 norm, if $^ 6 14 is a 
finite element solution of Poisson equation with Dirichlet boundary condition, we can find 
^h £ Wh so that <3>h + ^h is the finite element solution in Vh/2- In addition, suppose that (ph 
is the most accurate solution that a computer can compute in the sense that it would be out 
of memory when computing a finite element solution <&h/2 in Vh/2 directly Since the size of 
the linear system associated with ^ is smaller than $>h/2, if the computer can solve ^h, we 
can add ^h to $/,, to get $>h/2 achieving the next level of accuracy. In this paper, we discuss 
how to compute ^h- We shall construct compactly supported basis functions and a global 
supported basis function i^h,k, k = 1, • • • ,N h which span Wh- V\fc' s are called prewavelets 
and is a linear combination of these iph,k s an d hence is called a prewavelet subsolution. 

Prewavelets have been studied for more than 10 years (cf. [9], [5]). There are many 
methods available to construct compactly supported prewavelets over 2D domains under 
the L 2 norm. That is Wh = Vh/2 Q Vh under L 2 norm, e.g., in a series of papers [6], [7], 
[8], [11], and [1]. In 1997, Bastin and Laubin ([2]) explained how to construct compactly 
supported orthonormal wavelets in Sobolev space in the univariate setting. See also [1] for 
biorthogonal wavelets in Sobolev space. In [14], Lorentz and Oswald showed that there is 
no compactly supported prewavelets in Sobolev space or under H 1 norm based on integer 
translations of a box spline over R 2 . Since continuous piecewise linear finite element can 
be expressed by using box spline -Bin, the result in [U] ruins a hope to find compactly 
supported prewavelets under H 1 norm. But this is not an end of story. It is possible to 
construct compactly supported prewavelets in a semi- norm in the univariate setting in [TO] . 
It is also possible to construct compactly supported prewavelets in H r norm over each nested 
subspace, but the union of these prewavelets over all levels fails to be a stable basis for a 
Sobolev space (cf. [12]). Our new question is if we can find a prewavelet basis with as few as 
possible global supported prewavelet functions. Our anwser is affirmative. That is, there is 
a prewavelet basis for Wh with only one global supported basis function under the H 1 norm 
over rectangular domains. Also it is possible to find a compactly supported prewavelet basis 
for Wh under the H 1 norm for Poisson equation over a triangular domain (cf. |13j). 

The paper is organized as follows: We first explain that the Dirichlet boundary value 
problem of Poisson equation can be converted into a Poisson equation with zero boundary 
condition. An explicit conversion will be given. Thus the H l norm is now equivalent to the 
Hi semi-norm. Then we introduce some notation to explain the weak solution of Poisson 
equation and its approximation to the exact solution. These explanations are well-known and 
given in the Preliminary section §2. In §3, we explain how to construct compactly supported 
prewavelets under semi-norm. In §4, we explain how to implement our prewavelet method 
for numerical solution of Poisson equation. Finally in §5 we present some numerical results. 
Our numerical experiment show that the time for computing a finite element solution by our 
prewavelet method is about half of the time by the standard finite element method using the 
direct method for inverting the linear systems. If using the conjugate gradient method for 
the linear systems for the finite element method, the prewavelet method is still faster than 
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for sufficiently accurate iterative solutions. 



2 Preliminary 

Let us start with a square domain £7 = (0, 1) x (0, f) G R 2 . Consider the Dirichlet boundary 
value problem for Poisson equation: 



-Au(x,y) = g(x,y), 
u{x,y) = h{x), for 
u(x,y) = f 2 (x), for 
u(x,y) = f 3 (y), for 
u(x,y) = f 4 (y), for 



(x,y) G 
y = and 
y = 1 and 
x = and 
x = 1 and 



< x < 1 

< x < 1 

< y < 1 

< y < 1 



Without lose of generality, we may assume that /i(l) = /2(1) = /3(f) = ^4(1) — /i(0) = 
/ 2 (0) = / 3 (0) = / 4 (0) = 0. Otherwise, letting A(0) = / 3 (0) = a u / 3 (1) = / 2 (0) = a 2 , 
/2(f) = /4(f) = 03, /4(0) = /i(l) = a 4 , we define h(x,y) = a 1 + (a 4 - a x )x + (a 2 - a x )y + 
(a 3 + a x — a 4 — a 2 )xy, and v(x,y) = u(x,y) — h(x,y). Then the above Dirichlet problem 
becomes to: 

-Av(x,y) =g(x,y), (x,y)ett 

v(x, y) — fi(x) — h(x, 0), for y = and < x < 1 
v(x,y) = f'2{x) — h(x, 1), for y = 1 and < x < 1 
v(x, y) = f 3 (y) - h(0, y), for x = and < y < 1 
v(x,y) = fi{y) - h(l,y), for x = l and 0<y<l 

which satisfy the above assumption. 

Now let w(x) = v(x,y)-x(U(y) - h(l,y)) - (1 - x)(f 3 (y) - h(0, y)) - y(f 2 (x) - 
h(x, 1)) — (1 — y)(fi(x) — h(x, 0)). Then w(x) satisfies the equation 



-Aw(x,y) = gi(x,y), 
w(x,y) = 0, 



(x,y) e 
(x, y) G dVt 



with 9l (x, y) = g(x, y) + ^[-x(f 4 (y) - h(l, y)) - (1 - x)(f 3 (y) - h(0, y))} + £,[-y(f 2 (x) - 

h(x,i))-(i-y)(h(x)-h(x,o))\. 

If we can find solution for w, it is easy to get u(x,y). In the remaining paper, we 
only consider the Poisson equation with zero boundary condition: 
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-Au(x,y) = g(x,y), 
u(x,y) = 0, 



(x,y) G Q 
(x, y) G <9f2. 



(1) 



Next we define 



Hq(H) = {v G L 2 {Vl) : (v, v) s < 00 and v(x, y) = 0, (x, y) G dVt}, 
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where the inner product (u,v) s is defined by 

i(x, y) dv[ 

dx dx dy dy 



1 f 1 du(x,y) dv(x,y) | du(x,y) dv(x,y) 
u i v )s— / / ^ r\ ~r~ a o. axay. 



o ./o 



By using Poincare's inequality, ||u|| s = a/ (u, u) s is a standard Sobolev norm for Hq(£1). 
Suppose u, v G i?o(^)- Integration by parts yields 



r-2 z-2 

= / / g(x,y)v(x,y)dxdy 



./0 
2 /-2 



JO 
2 z-2 



o Jo 



—Au(x, y)v(x, y)dxdy 

du(x,y) dv(x,y) du(x,y) dv(x,y) 

a a ' a a dxdy 

dx dx dy dy 



= (u,v) s . 

Thus, a weak solution u to (pQ) is characterized by finding u G Hq(Q) such that 

(u,v).= (g,v), WveH^Q). (2) 

The following result is well-known. For convenience, we present a short proof. 
Theorem 2.1. Suppose g G C(fl) and w G C 2 (fi) satisfy (d)). TTien ii zs weaA; solution of 

Proof. Let v G .f/g(fi). Then integration by parts gives 

(g,v) = (u,v) s 

rl du(x,y) dv(x,y) du(x,y) dv(x,y) 

I ~ ~~ CLObcLy 



./0 



dx dx dy dy 

i ,i 

/ —Au(x,y)v(x,y)dxdy 
o </o 

= (-Au(x,y),v). 

It follows that (g — (—Au(x,y)),v) = for all v G if^fi). That is, g = —Am and hence, u 
satisfies ([I]). □ 

Next we introduce continuous linear spline space on Q = [0, 1] x [0, 1]. For conve- 
nience, let Nj = (2 J ' — l) 2 and j > 1. Denote Xji = i = y^ for i = 1, ..,2 J ' — 1. Clearly, 
the lines segment of x = x 3 j and u = yjp. divide the square Q into Nj sub-squares. The 
diagonal going from down-left to up-right of each sub-square divides the sub-square into two 
congruent triangle. We will refer to the set of all such triangles as a Type-1 triangulation of 
fl (see Figure 1). 
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Figure 1. Type-I triangulation with j=2. 



Define <$ ik to be linear spline with support on the hexagon with following vertices 

%•(*._!)), (xji,yj( k -i)), (xj(i+i),Vj(k)), (xj(i+i),yj(k+i)), ( x j(i+i),yj(k)), ( x j(i-i),yj(k))} 

and fakiXji' ,y jk ') = 8i,i'8k,k'> where <5 M / = if i' ^ i and 1 if i' = i. 

Let Vj = span{(p J ik , i = 1, .., 2 J ' — 1, k = 1, .., 2 J — 1} be the subspace of Hq (Q). By 
following lemma, there exists a unique Uj £ Vj satisfying 

(u j ,v) s = (f,v) WveVj. (3) 

Uj is the standard finite element solution in Vj. The following result is well-known. For 
completeness, we include a short proof. 

Lemma 2.1. Given g £ L 2 (f2) ; /ias a unique solution. 

Proof. Reorder the basis functions (j)^ to <j) m , m = l,...,Nj and let Uj = X] a m < / , m- Denote 

k mn = (0m,0n)s and F m = (f,<fi m ) for m = 1, , iVj. Set A = (a m ) to be the coefficient 

vector, K = [fc m „]i< m) „<Ar to be the stiff matrix, and F = (F m ) to be the right hand side 
vector. Then the solutions in ([3]) is written in the following matrix equation form 

KA = F. (4) 

We claim that the solution for above equation always exists and is unique. Otherwise there is 

a nonzero vector c such that Kc = 0. Write c = (c m , m = 1, , Nj) and let v = X^=i c «0« 

be the linear spline. Then Kc = is equivalent to 

(v,(f) m )s = Vm = 1, • • • , Nj. 

Multiplying (v, (p m ) s by c m and summing over m yields (v, v) s — 0. Thus, v = a + bx + cdxy. 
Boundary condition implies v = 0. Since {4> m } are linear independent, c = and hence, the 
solution is unique. □ 
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Let us discuss the error between u and Uj. It is standard in finite element analysis 
(cf. |3J). For completeness we present a simple derivation. Subtracting ([3]) from ([2]) implies 



(u — Uj, w) t 







\/w G Vj. 



(5) 



Then for any v G Vj 



\U — Uj 



(u — Uj, U — Uj) s 

(u — Uj, U — v) s + (u — Uj, V — Uj) t 
(U - Uj,U- V) s 



< \\U — Uj\\ s \\U — v\\ s 



It follows that 1 1 it — Uj\\ s < \\u — v\\ s for any v G Vj. Thus we have proved the following. 
Lemma 2.2. ( Cea's Lemma) \\u — Uj\\ s = min{||-u — v\\ s : v 6 Vj}. 

Given u G C°(fl), let Uj G Vj be the interpolant of v : 



The following error estimate is well-known. 
Lemma 2.3. Suppose u G C 2 (f2). T/ien 



\U — Uj\\ s < 



12 



2-i 



d 2 u 



dx 2 



du du 
dx dy 



d 2 u 



dy 2 



Proof. The proof is elementary and is left to the reader. See [13] for detail. 



□ 



3 Multiresolution and Prewavelets over Type-I trian- 
gulations 

We start with the definition of multi-resolution approximation of Hq(Q): 

Definition 3.1. A multiresolution approximation of Hq(Q) is a sequence of finite dimensions 
subspaces Vj, j G Z + of Hq(Q) such that 

(1) VjCV j+1 , jeZ+ ; 

(2) U°l i Vj is dense m H*(Q). 
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Let T- 7 be the type-1 triangulation with 2Nj triangles. Naturally, let T^ 1 be the 
uniform refinement of P. Let Vj be the continuous piecewise linear spline space defined on 
the previous section. That is, Vj = span{(f> J ik , i = 1, .., 2 J 1 — 1, k = 1, .., 2 J — 1}, where (f> J ik are 
continuous piecewise linear functions which is 1 at (xji,yjk) and zero at all other vertices. 
Let Vj+i = span^i^ 1 ,i = 1, .., 2 J+1 — 1, k = 1, .., 2 J+1 — 1}, and (xj + i^, yj+i,k) are the vertices 
on the j + 1 level Type-1 triangulation. Then the refinement equation is easily seen to be 

jj +-6 j+1 +-6 j+1 +-6 j+1 +-6 j+1 +-6 j+1 +-6 j+1 

Vik ~ ^2i,2k ' 2^ 2i - 1 ' 2k 2 2i ~ 1 > 2k ~ 1 2 2i ' 2k ~ 1 2 2i+1 ' 2k 2 2i+1 ' 2k+1 2 2i ' 2k+1 ' 

See the Figure 2. 




(0,0) (1/2^,0) 

Figure 2. Dilation relations 

The main purpose of this paper is to build a basis for the orthogonal complement Wj 
of Vj in Vj+i under the inner product (•, -) s . Suppose we have the Wj. Then Vj+i = Vj + Wj 
under the Hq{VL) inner product. For a solution Uj satisfying (3), we do not have to find out 
the solution for 

Uj + i G Vj+i such that (v,j + i,v) s = (g,v) Vf G Vj + \. 

Instead, we only need to find solutions for 

Wj G Wj such that (wj,v) s = (g,v) Vi> G Wj. 

Then we have Wj + Uj = Uj +1 . Ideally, we hope the supports of basis functions for Wj are 
small, since small support can accelerate the calculations of (g,v) s . As explained in the 
Introduction, there is no compactly supported prewavelets for Wj. Neverthless, we shall 
construct basis functions with only one globally supported basis function for Wj in the 
following. 

Clearly the Tj can be continuously refined and hence we will have a nested sequence 
of subspaces 

V 1 C V 2 C V 3 C V 4 C V 5 

to span Hq(Q) by Lemmma 2.3 since C 2 (Q) is dense in Hq(Q). 
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Let Wj C Vj + % be the orthogonal complement of Vj in Vj+x for each refinement level 

j i.e. ; 

Then we get the decomposition 

for any j > 1. The weak solution to the Poisson equation ((T]) at Vj+i can be built by 

ttj+1 = Ui + Wi + u> 2 H 1" Wj- 

We now focus on building basis functions for the orthogonal complement Wj. By 
direct calculation, we obtain the following lemma immediately 

Lemma 3.1. We have (<^ fc , <frg* h , ) s = 2, 

(4. 4t-i,2k, )s = 1/2, (4, 0iS*-i, > s = 1/2, (4, 0^ lj2fc , ) s = 1/2, 

(0ijfc, 02i+l> >- = V 2 ; (0ifc, 0itl,2fc-l, )* = ^ (0* <%i+l,2k+l> )« = ^ 

(4, >* = -1/2, (04, 0L + + 2 , 2 *, >, = -1/2, (4> 0£i-2, >* = -1/2, 

0^+2' = "I/ 2 , (0ifc, 02^-2,2^-2, )- = °. (0ifc, <Pit+2,2k + 2^ )s = 0, 

(0i, 02^-2,2^-1. = -I/ 2 , (0L 02^-1,2^+1, >* = - 1 ; (0L ^«+i l 2fc+a» )* = -I/ 2 , 

(0ijfc, ^2,2^+1, = -I/ 2 , (0ifc, 0S- 1 l,2fc-l, >- = "I. (0ifc, 02^-1,2^-2. ) S = -I/ 2 , 

(0ifc; 0i' + fc 1 '; )s = 0; f or other i , k' which are not listed above. 

Let ifi* be a function in Wj. Since Wj C Vj+i, let us write ipi = ^2 ik b^ for some 
unknown coefficients b^. Then by orthogonal condition (<f>l, k f, t/> j ') s = 0, we need to solve the 
following equations. 

Each (i A;') determines one equation. Since there are Nj elements in the set Vj, they determine 
the Nj equations. These Nj equations with Nj+\ coefficients, b^. There are at least Nj + i—Nj 
degrees of freedom. The solution space of these equation system should be the Wj. The 
linear independence of 0^, k , implies that the coefficient matrix of the above linear system 
is of full rank. Hence, there are Nj + \ — Nj linear independent solutions which constitute a 
basis for Wj. 

Definition 3.2. Let V J ™ 1 = span{(/P ik , 2 = 1,.., 2m — 1, k = 1, .., 2m — 1} be a subspace of 
Vj+i. Let WJ 1 be subspace ofWj such that W™ = Wj^V^. 

Obviously C V} +1 C V? +1 C . . . C V? +1 = V j+1 , and Ic^cHjc.C Wf = 
Wj. There is no nonzero solution of (Q in space of V^i- However, there are five solution of 
(JBJ) in space V? +1 . They are solutions of the following system of linear equations. 
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E b * <4 +1 ^ii> s = o, E 

l<i,fc<3 l<i,fc<3 
l<i,fc<3 l<i,fc<3 

They are equivalent to the following equations. 



tit)* 

4f)s 



^1,2, 



^2, 



0. 

0. 



/ Ki \ 

&2,1 
&3,1 
&1,2 
&2,2 
&3,2 
&1,3 
&2,3 
V fo 3,3 / 











v o y 



Using Lemma 13.11 we obtain the following equations. 



/ 1 


1/2 


-1 


1/2 


2 


1/2 


-1 


1/2 


1 \ 







-1/2 


1 





-1/2 


1/2 








-1 



















-1/2 





-1/2 


1 












-1/2 


-1/2 





1 


1/2 


- 1 1 





( h A \ 

h,i 
h,2 
b 2 ,2 
h,2 
h,3 
h,3 

V 6 3,3 / 











v o y 



The rank of the left matrix is four, because 



dent. So there are five solutions as shown below. 

/ 



^2,2i 



are linear indepen- 



/ Ki \ 

b 2 ,i 
h,i 
h, 2 
b 2 ,2 
h, 2 
bi, 3 
b 2 ,3 

V & 3,3 j 





2 


1 


V o / 



or 



/0\ 

2 
1 






v o y 



or 



1 
1 

1 

-1 




V o / 



or 








\ 




/ 





\ 













-l 































1 






-1 




or 











1 








-1 






















1 








1 






1 


/ 




V 





/ 
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More precisely, 



= + as shown in Figure 3; (7) 

ift* = 2^ + as shown in Figure 4; (8) 

^i;? = -02? + 03? + 02? + 03? as shown in Figure 5; (9) 

^'j = + + - 0g x as shown in Figure 6; (10) 

= 0i? + 02? — 02? — 03? as shown in Figure 7. (11) 




(0,0) (1/2', 0) (2/2'', 0) (0,0) (1/2', 0) (2/2', 0) 

Figure 3. Figure 4. 




(0,0) (1/2', 0) (2/2'", 0) (0,0) (1/2', 0) (2/2', 0) 

Figure 5. Figure 6. 
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(0,2/2*) 



(0, 1/2*) 

















1 / 





(0,0) 



(1/2*, 0) (2/2*, 0) 
Figure 7. 



Now we consider Vf . Similarly, there are 25 non-zero coefficient for linear system ([6l) 
and the coefficient matrix of rank 9. So the dimension of solution space of Wj is 25 — 9 = 16. 
The first five of them are the same to the wavelet functions in (7)- ffTTl) . The other 11 are 
given below. 



-0; 

l/j: 
ip: 
if;: 
if): 



+ 



4? 

6 j+1 

^1,4 

0S 1 



& + 



6* +1 

9 3,4 
b j+1 

6* +1 

9 4,3 



^4,3 



6 j+1 

^2,4 
^4,4 

6 j+1 

^4,2 
^* +1 

r4,2 

6 j+1 
r4,4 

^* +1 
r2,4 

^3,4 



as shown in 
as shown in 
as shown in 
as shown in 
as shown in 
as shown in 
as shown in 
as shown in 
as shown in 
as shown in 
as shown in 



Figure 8; 
Figure 9; 
Figure 10 
Figure 11 
Figure 12 
Figure 13 
Figure 14 
Figure 15 
Figure 16 
Figure 17 
Figure 18 




(0,3/2*) 
(0,2/2*) 
(0, 1/2*) 



































/ 


2 1 

/ 







(0,0) (1/2*, 0) (3/2*, 0) 



(0,0) (1/2*, 0) (3/2*, 0) 



Figure 8. 



Figure 9. 
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(0,3/2'") 
(0,2/2*) 
(0, 1/2') 















L 1 










1 1-7 



















(0,0) (1/2', 0) (3/2', 0) 
Figure 10. 



(0,3/2') 
(0,2/2') 
(0, 1/2') 



























/ 

L 1 










1 \y 







(0,0) (1/2', 0) (3/2', 0) 
Figure 12. 



(0,3/2') 
(0,2/2') 
(0, 1/2') 















^1 - 










1 


1 —7 

















(0,0) (1/2', 0) (3/2', 0) 
Figure 14. 



(0,3/2') 
(0,2/2') 
(0, 1/2') 

















L 1 










1 ly 

















(0,0) (1/2', 0) (3/2', 0) 
Figure 11. 



(0,3/2') 
(0,2/2') 
(0, 1/2') 

























/ 1 


1 








/ \/ 







(0,0) (1/2', 0) (3/2', 0) 
Figure 13. 

(0, 3/2') XXX/ 
(0,2/2')/^ / 
(0, 1/2') /l Y / / 

(0,0) (1/2', 0) (3/2', 0) 
Figure 15. 
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(0,3/2'") 
(0,2/2*) 
(0, 1/2') 























l / 




















(0,0) (1/2', 0) (3/2', 0) 
Figure 16. 



(0,3/2') 
(0,2/2') 
(0, 1/2') 



(0,0) (1/2', 0) (3/2', 0) 
Figure 17. 



































/ 7 
/ 1 / 







(0,0) (1/2', 0) (3/2', 0) 
Figure 18. 

The above computation can be carried out on VJ 1 for n = 3, ....2' — 1. We have thus 
obtained five types of wavelet functions: 

is supported next to the vertical boundary and is called vertical boundary wavelet. 

Co = + 

called horizontal boundary wavelet, is supported next to the horizontal boundary. The next 
three types are supported inside the domain. The following 

Vi,k ^i+l,fe+l ^ Vi + 2,fc+l ^ S t i + l,fc+2 ^ V^ + 2,fc+2 

is called interior wavelet of first kind. We call 

ib jA - -fh j+1 + rh j+1 4- r/>' +1 4- rh j+1 
interior wavelet of second kind. The last one 



j'+i 

2i-l,2fc 



+ 



^2i,2fc+l 



J 1i,1k-\ 



J 2i+l,2k 



is called interior wavelet of third kind. 
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Theorem 3.1. All the five types of wavelets in the are linear independent for 1 < n < 
2j — 1. That is, for each 1 < n < 2j — 1, the following functions 

V'o'i k = l,..,n-l, 
k = l,..,n-l, 
1 < i, k < n - 1, 
l<i,fc<n-l, 
l<i,k<n-l 

are linear independent. 

Proof. Let us prove it by induction. It is true for n = 2 and for n — 3. Suppose it is true 
for n = p, that is, 





fc = 


-- l,..,p- 


-i; 




k = 


-- l,..,p- 


- 1; 


¥%fc> 


1 < 


i,k < p 


- 1 




1 < 


i,k < p 


- 1 




1 < 


i,k < p 


- 1 



are linear independent. For n — p + 1, there are 6p — 1 new functions which are 

V'ot fc = PI 
^k,o> ^ = PI 

^(fe, i or k = p; 

V'i'fe, i or fc = p; 

V'-J i or fc = p. 

Suppose they are not linear independent. That is, one can find 



aj fc , i or k=p; 

af fe , i or fc = p 

such that 

«vs+a 2 Co+ e "'''./■'• ''./.•' • e <^a+ e <^a+^'=o, (12) 

i or fc=p i or k=p i or fc=p 

where tp' is linear combination of the following functions: 



Y0,ki 


fc = 


~- 1, ..,p~ 


i; 


lb j ' 2 


fc = 




i; 




1 < 


i,k < p — 


l 




1 < 


i,k < p — 


l 




1 < 


i,k < p — 


l 
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By the definition, (fyt+i 2k+n * = P or k = p appear only once in ^'jjj, i = p or k = p , iffi and 
V'p o- Since are linear independent, that is, a? fc = 0, % or k = p, a 1 = 0, and a 2 = 0. 
Thus the equation ([T2]) can be simplified to 

E + E <>■■'■■'! • °- ( 13 ) 

i or fc=p i or k=p 

By the similar reason, 4> J 2 ^l k , i = p or k = p appear only once in ipf' k , i = p or k = p. Since 

(ftik 1 are linear independent, a\ k = 0,iork = p. Thus the equation ffT3l) can be further 
simplified to the following equation 

E «W3+v>' = o. 

i or A;=p 

Similarly, af A = 0,i or k = p too. Thus the equation (TT2]) is reduced to 

V>' = 0. 

By induction hypothesis, all the coefficient of ip = are zeros. Hence, 

V'qt fe=l,..,n-l, 
V'feJ fc=l,..,n-l, 
^(f, 1 < i, A: < n - 1, 
1 < i, fc < n - 1, 
V'i.f > l<i,A;<n — 1 

are linear independent. □ 

Theorem 3.2. AZ/ i/ie /iue types of wavelets in the W™ form a basis of WJ 1 for 1 < n < 
2 j — 1 . That is, 

W? = span{^ Vg, V#, 1$, 1$, 1 < *, * < * " U 

for 1 < n < 2j — 1. 

Proof. The dimension of Wj 1 is (2n — l) 2 — (n) 2 = 3n 2 — 4n + 1. It is easy to count that 
there are (2n — l) 2 — (n) 2 = 3n 2 — 4n + 1 functions in the following set 

^JS, 1 < h k < n; 

1 < /c < n; 
V'i'fc, 1 < A; < n 

which all belong to the space W™. Since they are linear independent, they form a basis for 
space Wj 1 , where 1 < n < 2j — 1. □ 

15 



Finally we need to find wavelets in Wj J \Wj J ~ l . The computations are the same to 
the above except for that there is one globally supported basis function. In fact the following 
pictures show the basis functions located on the top boundary of the domain Q. (We omit 
the pictures for the basis functions on the right vertical boundary which are symmetric with 
respect to the line y=x are those basic functions on the top horizontal boundary of Q.) 




Figure 23. 



16 



Figure 24. 



The last one (cf. Figure 24) is the only special basis function since it is not local supported. 
The numbers of all these wavelets in Wj^Wj 3 " 1 amount to 2 J+3 — 8 which is equal to the 
number of dimension of V^-^V^ 1 . 

Theorem 3.3. All the wavelets in the W^XWj 3 ' 1 are linear independent and form a basis 
for V^l^V^ 1 which is spanned by the functions in {0f£\ 2 J+1 — 2 < i, k < 2- ?+1 — 1}. 

Proof. Let us just concentrate on the basis functions in V^-^V^ 1 and in Wj 7 \Wf ' _1 ■ Then 
the scaling matrix between two sets of basis functions is the following matrix up to a constant 



/ D 

Bl B2 

Bl B2 

Bl B2 



A 



C3 C3 C3 



Bl B2 

CI C2 
C3 C3 



C3 
C4 

B2' BY 

B2' BY 



\ 



B2' BY 
D> J 



where 



D = ( 1 2 ) , Bl 



( 1 2 \ 

1 -1 

1 1 

V 1 - 1 / 



B2 = 



( \ 


0-100 

\1 1 0/ 
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D' = ( 2 1 ) , BY 



/-ll \ 

1 1 

-1 1 

V 2 1/ 

/ \ 
10 
0-1-10 

V i i oo/ 

C3 = ( 1 ) . 
Let E = (m n 0). By the row operations we have 

/ m n 



52' = 



/ 1 1 \ 

0-10 



\ / 



CI = 



( 1 2 \ 

1 -1 

1 1 

V 1 - 1 / 



C2 



CA = ( 2 1 ) , 




2 

1 
1 1 

1 



0-1 
-1 1 1 
1 



2 

1 



1 



/ 



m n 



1 10 0-1 
1-11 1 / 



-n 2m 

2m —n 
n 



m 

2m + n 2n 
1 
1 
1 



V 

Similar for B' . Thus by row operations, 

^3 G 3 





2 

-1 

1 
1 -1 1 



-1 
1 J 



A 



A 



21-2 



G 2 j-2 



C' 



G' 



23-2 ^2i-2 



A' 



G' 2 A> 2 



G[ A[J 
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where A n is an upper triangular matrix of size 4x4 while A' n is a lower triangular matrix 
of size 4x4 which are given below. 



A 1 



( 1 


2 












-1 


1 











1 


-1 


V 








2 




(o 








o\ 
































V i 








0/ 



/ \ 





\ i o o o / 



A, 



( 1 1 \ 

-12 

2 -1 

V 1 J 



A' 



i A r 



( 2 

-1 n 

n 

\ 



/ n 
\ 
1 



2 \ 

-1 71 

n -1 
2 / 



, C ra — 



\ 



G' 



2 n / 

and the matrix (C[ C 2 ) is the following matrix 



/ n \ 





\ J 



( 2? 



+i 



2 

1 
1 



2^' +1 - 5 



\ 



-1 


-1 







1 
1 



1 



2 



V 



o 

2 



2i +l - 5 / 

It is easy to see the rank of (C[ C 1 ^) is 8. Thus the rank of A is 8(2 J ) 



/ \ 





\ n / 



prewavelet functions constructed above in the Wf 3 \Wf J ' 



form a basis of V^ 2 | ^ 1 \V J+1 



2^-1 



3. Thus, all the 
1 are linear independent and hence 

□ 



It is easy to see that the coefficients of the prewavelet functions in W, 



2J -1 



in terms 

2P-1 



of the basis functions of V^XV^ 1 are all zeros. Thus the prewavelet functions in Wj 

together with the prewavelet functions in V^XV^ 1 are linear independent. It follows the 
main result in this paper. 

Theorem 3.4. All the prewavelet functions in the Wj^W^' 1 and the prewavelet functions 
in Wf~ x form a basis for Wj. 
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4 The Prewavelet Method for Poisson Equation 



Let us use the basis functions of Vj and Wj to solve Poisson equation ([T]). Mainly we explain 
how to compute hj G Wj. Let gj G Vj and G Vj + \ be two FEM solutions. We aim to 
show that hj + gj = gj+i- 

By a reordering the indices (i,k),l < i,k < 2- ? in a linear fashion, let Vj = 
span{(j){, , (f) 3 N .}- Also, we reorder all five type wavelet functions as well as the glob- 



ally supported wavelet to denote Wj 
vectors, 



span{if){, ip 



N 



3 + 1 



-N. 



}. Let W be following 



\ <t>k 1 

Then we have the following equations 



3+1 



-JV, 



+1 



where i?j is Nj x iVj+i refinable matrix, and Cj is a wavelet matrix of size (Nj+i — Nj) x Nj + \. 
Let and i?j be the following matrices: 



(4 Ah 



J 2i VAT 



E; 



( 



MM)* 



A' 



3+1' 



-AT,- 5 



A' 



3+1' 



-AT.) 



N i+X -Njl» 



it 



It is easy to see that BjDj + iCj = is equivalent to Vj-LWj. Clearly, we have Dj = 
BjD j+1 Bj and Ej = CjD j+1 Cj. 

Let gj be the projection of g in Vj, and hj be the projection of g in Wj. Since 
Vj Wj = Vj+i, gj+hj will be equal to gj+i- Let us write gj = J2j=i a «0; = ( a i> a 2> ■• ■> OiV,)^- 



Similarly, /ij 



t,6 2> ....,6 



iv 



3 + 1 



-AT, 



l^/- 7 , and gj + i = (ci,c 2 , CAr. +1 )$ J+1 . By computing the 



weak solutions hj,gj, and in Wj-, Vj, and VJ+i, respectively, we have 



( d\ \ 

V a ^ ) 



( 



\ 
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It follows 



( 61 \ 

b 2 



c 2 



\ (^, +1 -AT,^> / 



/ <0i + \<?> \ 
(0 J 2 +1 ,^) 



/ ai \ 

a 2 



V a ^ J 



( Cl \ 

C 2 

V Ctv j+1 J 



(Ej) 



-1 



/ (4,9) \ 
(4,9) 



\ (4 Nj ,9) ) 

( (ti,9) ■ 
(4,9) 

\ (^N j+1 - Nj ,9) J 



Vj + l 



/ (4 +1 ,g) \ 



(D^r 1 



The above linear systems provide a computational method to find gj, hj. 

We now show hj + gj = gj+i- That is, gj + i can be computed by using hj and gj 
only. Indeed, we have 

9j = (ai,a 2 , ,a Nj )®> = ($ J ) T (ai,a 2 , ,a nj ) T 

= (& +1 ) T Bj( ai ,a 2 , ,a N f 

= (& + yBjD-\(4,g), (4,g), • • ■ (^,g)f 



3 3 



((& +1 )) T Bj(B j D j+1 Bjy 1 B j ((^ 1 ,9), (4 +l ,g), ■ ■ ■ (4J'g)) 



3+1 ■ 



and 



Similarly, 



hj = ((^fcjiCjD^cjy'c^^g), (4 +1 ,g), ■ ■ ■ (4j'g)) 



7+1 



j'+i' 
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In order to show hj + gj = gj+i, we only need to prove 

Bj(B j D j+1 Bjy 1 B j + C< [C,!),. : C<) V' = D£ v (14) 

Notice that Bj and Cj are not square matrices. That is we can not invert Bj and Cj. 
Consider 

- V C,D J+1 CJ J 

by using the orthogonal conditions of Vj and Wj. Then we have the following equation 

V CjD j+1 J ( 1 * j ^ o (CjD^cjy 1 

where / stands for the identity matrix. In other words, we have 

( IfDZ ) ( mW^XfT 1 CJ(C 3 D J+1 CJ)- 1 ) = / 
which can be rewritten in the following form 

( Bj(B jDj+1 Bj)-i Cj(CjD j+1 Cj^ ) ( B c f D ^ ) = I. 

Hence we have 

Bj(BjD j+1 Bj) BjD j+1 + Cj(CjD j+1 Cj) CjD j+1 = I 

or 



BjiBjD^Bjy'Bj + cjiCjD^cjy'cj = n.y. 

which is (fl4"]l and hence hj + = gj+i- 



5 Numerical Experiments 

We have implemented the prewavelet method for numerical solution of Poisson equations 
over rectangular domains in MATLAB. We would like to demonstrate that our prewavelet 
method is more efficient than the standard FEM method. 

In the following we provide three tables of CPU times for numerical solutions based 
on our prewavelet method and the standard finite element method for various levels of 
refinement of an initial triangulation (Tq which consists of two triangles) of the standard 
domain [0, 1] x [0, 1]. 
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Let Vj be the continuous linear finite element space over triangulation Tj which is 
the jth refinement of T . For a test function u which is the exact solution of Poisson equation 
([1]), the finite element method is to compute Uj G Vj directly while our prewavelet method 
computes Uj by computing Wk, k = 1, • ■ ■ , j, i.e., Uj — U\ + W\ H h Wj-\. 

In the following we present three tables of CPU times for computing numerical 
solutions Uj,j = 4,5,6 for three test solutions by using these two methods. Note that we 
use the direct method coded in MATLAB to solve the associated linear equations. We 
shall present tables of CPU times based on Conjugate Gradient Method for the systems of 
equations next. 

For an exact solution u(x,y) = sin(27rx) sin(27n/) which clearly satisfies the zero 
boundary conditions, we list CPU times for computing numerical solutions Uj, j = 4, 5, 6 by 
using these two methods in Table 1. 

Table 1. CPU times to compute Uj by the two methods 





FEM method 


Prewavelet Method 


j=4 


0.164531 seconds 


0.204067 seconds 


j=5 


0.593587 seconds 


0.519293 seconds 


j=6 


13.960323 seconds 


6.222679 seconds 



For an exact solution u(x,y) = xy(l — x)(l — y), the CPU times for numerical 
solutions by these two methods are given in Table 2. 

Table 2. CPU times for computing Uj by the two methods 



CPU time 


FEM method 


Prewavelet Method 


j=4 


0.150836 seconds 


0.218282 seconds 


j=5 


0.574085 seconds 


0.558071 seconds 


j=6 


13.896825 seconds 


6.202557 seconds 



We list the CPU times for computing numerical solutions Uj, j = 4, 5, 6 of u(x, y) = 
xy(l — x)(l — y)e 8xy by using these two methods in Table 3. 

Table 3. CPU times for computing Uj by the two methods 



CPU time 


FEM method 


Prewavelet Method 


j=4 


0.144159 seconds 


0.186389 seconds 


j=5 


0.584828 seconds 


0.459181 seconds 


j=6 


13.877403 seconds 


6.139101 seconds 



It is clear from these three tables that the prewavelet method is much more efficient. 
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Next we use the Conjugate Gradient Method to solve the linear systems associated 
with FEM. Let us consider iterative solution to Uj for j — 6 with various accuracy. First let 
us consider the exact solution u(x,y) = sin(27ra;) sin(27rc/). 

Table 4. CPU times for approximating the FEM solution w 6 by Conjugate Gradient Method 



e 


CPU times 


1(T 8 


5.411852 seconds 


10- 9 


5.783497 seconds 


10 -io 


6.221683 seconds 


10- 11 


6.616816 seconds 


10- 12 


6.917468 seconds 


10- 13 


7.836775 seconds 



To approximate the FEM solution u e of the exact solution u(x, y) = xy(l — x)(l — y) 
by the Conjugate Gradient Method, we list the CPU times in Table 5. 

Table 5. CPU times for approximating the FEM solution Uj by Conjugate Gradient Method 



e 


CPU times 


10- 8 


4.476794 seconds 


10- y 


4.878259 seconds 


10 -io 


5.306747 seconds 


10- 11 


5.887849 seconds 


10- 12 


6.811317 seconds 


10- 13 


6.754465 seconds 



Finally let us consider the CPU times to approximate the FEM solution uq of u(x,y) = 
xy(l — x)(l — y)e 8xy by the Conjugate Gradient Method. 

Table 6. CPU times for approximating the FEM solution by Conjugate Gradient Method 



e 


CPU times 


10- 8 


10.110517 seconds 


10- y 


10.740035 seconds 


10 -io 


11.319618 seconds 


10- 11 


11.810142 seconds 


10- 12 


12.320903 seconds 


10- 13 


13.103407 seconds 



It is clear from all six tables, if we want an accurate iterative solution of u§ within 10 
the prewavelet method appears better. 
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